﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using DMC=DiffusionMonteCarlo;

namespace Examples
{
    public class MorseOscillator
    {
        private string result;
        public string Result
        {
            get { return result; }
        }
        public void Compute()
        {
            DMC.Input input = new DMC.Input();

            input.NumberOfParticles = 1;
            input.SpaceDimension = 1;

            input.Seed = 22312345L;
            input.InitialNumberOfReplicas = 1000;
            input.RunTime = 1000.0;
            input.ValueOfTimeStep = 0.05;

            double[,] minmax = { { -20.0, 20.0 } };
            input.LimitsOfCoordinates =  minmax;
            double[,] ip = { {0.0} };
            input.InitialPosition = ip;
            double[] delx = { 0.1 };
            input.ValueOfLengthInterval = delx;
           

            DMC.PotentialEnergy<double[,], double> V
                = new DMC.PotentialEnergy<double[,], double>(Function);
    
            DMC.PathIntegral pi 
                = new DMC.PathIntegral(input, V);
            
            DMC.Output output = pi.Compute();

            result += output.GroundStateEnergy.ToString() + "\r\n";
            for(int i = 0; i < output.GroundStateWaveFunction.Length; i++)
                result += output.GroundStateWaveFunction[i].ToString() + "\t";

            double[,] xy = new double[3, output.GroundStateWaveFunction.Length];
            for (int i = 0; i < output.GroundStateWaveFunction.Length; i++)
            {
                xy[0, i] = input.LimitsOfCoordinates[0, 0]
                    + ((double)i + 0.5) * input.ValueOfLengthInterval[0];
                xy[1, i] = output.GroundStateEnergy;
                xy[2, i] = output.GroundStateWaveFunction[i];
            }
            Graph graph = new Graph(xy);
            graph.Show();
        }

        public double[,] xy;

        private double Function(double[,]  x)
        {
            return 0.5 * (Math.Exp(-2.0*x[0,0]) - 2.0*Math.Exp(-x[0,0]));
        }
    }
}